Introduction to R

본 강의의 목표

• R 프로그래밍을 시작하고 R 문법에 익숙해진다.
  
• 데이터를 입력받아 다룰 수 있다.
  
• 다양한 Table, Plot들을 만들 수 있음을 안다.

“ R is a free software environment for statistical computing and graphics. It compiles and runs on a wide variety of UNIX platforms, Windows and MacOS. ”

https://cran.r-project.org/bin/windows/base/ & https://www.rstudio.com/products/rstudio/download/

R 기초연산

벡터(vector)

R의 기본 연산 단위는 벡터이다.

x=c(1,2,3,4,5,6)            ## vector of variable
y=c(7,8,9,10,11,12)
x+y; x*y; sqrt(x)

## [1]  8 10 12 14 16 18

## [1]  7 16 27 40 55 72

## [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490

mean(x); sd(x); length(x)

## [1] 3.5

## [1] 1.870829

## [1] 6

x[2]; x[-2]; x[c(1,3,5)]

## [1] 2

## [1] 1 3 4 5 6

## [1] 1 3 5

벡터만들기

벡터를 생성하는 다양한 방법

## Sequence
v1=seq(-5,5,by=.2); v1 

##  [1] -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2
## [16] -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2  0.0  0.2  0.4  0.6  0.8
## [31]  1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4  2.6  2.8  3.0  3.2  3.4  3.6  3.8
## [46]  4.0  4.2  4.4  4.6  4.8  5.0

## Repeat 
v2=rep(1,3); v2

## [1] 1 1 1

v3=rep(c(1,2,3),2); v3              ## Repeat for vector

## [1] 1 2 3 1 2 3

v4=rep(c(1,2,3),each = 2); v4       ## Repeat for vector : each

## [1] 1 1 2 2 3 3

for, if, else, ifelse 문

1

## for loop
for (i in 1:3){
  print(i)
}

## [1] 1
## [1] 2
## [1] 3

i=0
for (j in c(1,2,4,5,6)){
  i=i+j
}
i

## [1] 18

2

## if 
x=5
if (x >=3 ){
  x=x+3
}
x

## [1] 8

x=5
if (x >=10){
  print("High")
} else if (x >=5){
  print("Medium")
} else {
  print("Low")
}                                ## if, else if 주의: 반드시 } 와 같은 줄에 위치하도록.

## [1] "Medium"

## ifelse
x=5
y=ifelse(x==5,"OK","Suck")       ## ifelse(조건,참일때,거짓일때)
y

## [1] "OK"

Apply 문 : apply, sapply, lapply

1

R은 벡터 기반의 연산을 지원하므로 이를 이용하면 여러가지 계산을 간단히 할 수 있다.
apply는 행 또는 열 단위의 연산을 쉽게 할 수 있도록 지원하는 함수이다.

mat=matrix(1:20,nrow=4,byrow=T)   ## 4행 5열, byrow=T : 행부터 채운다. 
mat

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    2    3    4    5
## [2,]    6    7    8    9   10
## [3,]   11   12   13   14   15
## [4,]   16   17   18   19   20

apply(mat,1,mean)                 ## 1: 행- 행별로 평균

## [1]  3  8 13 18

apply(mat,2,mean)                 ## 2: 열- 열별로 평균

## [1]  8.5  9.5 10.5 11.5 12.5

2

sapply는 실행결과를 벡터로 출력하고 lapply는 실행결과를 리스트로 출력한다.

sapply(1:nrow(mat),function(x){mean(mat[x,])})       

## [1]  3  8 13 18

sapply(1:ncol(mat),function(x){mean(mat[,x])})

## [1]  8.5  9.5 10.5 11.5 12.5

lapply(1:nrow(mat),function(x){mean(mat[x,])})

## [[1]]
## [1] 3
## 
## [[2]]
## [1] 8
## 
## [[3]]
## [1] 13
## 
## [[4]]
## [1] 18

unlist(lapply(1:nrow(mat),function(x){mean(mat[x,])}))     ## Same to sapply

## [1]  3  8 13 18

3

for문의 과도한 사용은 분석의 속도를 저하시키는 가장 큰 요인 중 하나이다.
따라서 for문을 최대한 덜 쓰고 벡터연산을 활용하는 것이 실행시간을 줄이는 핵심이다.

## for 문을 이용한 합구하는 함수 만들기 
sum_f=function(x){
  out=0
  for (i in 1:x){
    out=out+i
    }
  return(out)
}
system.time(sum_f(10^6))                  ## system.time - 시간 잰다.

##    user  system elapsed 
##   0.048   0.000   0.048

system.time(sum(as.numeric(1:10^6)))      ## 내장 sum 함수- 더할 것이 많으면 as.numeric을 넣어줘야 한다. 

##    user  system elapsed 
##       0       0       0

데이터 불러와서 작업하기

데이터 불러오기

데이터를 불러오기 전에 미리 디렉토리를 지정하면 그 다음부터는 편하게 읽고 쓸 수 있다.
여기서 주의해야 할 점은 폴더간의 구분을 / 로 해야 한다는 점이다. R 은 유닉스 기반의 언어이기 때문이다.

setwd("/home/heejooko/ShinyApps/introduction to R")
getwd()

## [1] "/home/heejooko/ShinyApps/introduction to R"

예제로 파일을 불러와 살펴보자.

a<-read.csv("practice.csv")
a

읽은 데이터 살펴보기

1

head(a); tail(a)

a$SEX <- ifelse(a$SEX=='0',NA,a$SEX)
tail(a)

2

class(a)

## [1] "data.frame"

names(a)

## [1] "CASEID"    "AGE"       "SEX"       "HEIGHT"    "WEIGHT"    "DM"       
## [7] "HTN"       "PH"        "TREATMENT"

nrow(a); ncol(a)

## [1] 50

## [1] 9

a[,1:5]

a[20,]

3

magrittr 패키지의 %>% 연산자를 이용하면 직관적인 코드를 작성할 수 있다.

# install.packages("magrittr")
library(magrittr)
a %>% head                ## head(a)와 같다.

a %>% head %>% nrow       ## nrow(head(a))와 같다.

## [1] 6

4

키에 대한 통계량을 살펴보자

mean(a$HEIGHT)                                 ## 결측치 없을 때만

## [1] NA

mean(a$HEIGHT, na.rm=T)                        ## 결측치 빼고

## [1] 162.6043

round(mean(a$HEIGHT, na.rm=T), 1)              ## 반올림

## [1] 162.6

a$HEIGHT %>% mean(.,na.rm=T) %>% round(.,1)    ## 위와 같은 명령

## [1] 162.6

5

아래의 함수 모두 na.rm=T을 빠트려서는 안된다.

sd(a$HEIGHT,na.rm=T)                          ## 표준편차

## [1] 18.66815

var(a$HEIGHT,na.rm=T)                         ## 분산

## [1] 348.5

median(a$HEIGHT,na.rm=T)                      ## 중간값

## [1] 163.9

IQR(a$HEIGHT,na.rm=T)                         ## 25%-75% range

## [1] 14.1

quantile(a$HEIGHT, na.rm=T)                   ## quantile

##    0%   25%   50%   75%  100% 
##  51.3 158.4 163.9 172.5 183.0

max(a$HEIGHT,na.rm=T)                         ## Max

## [1] 183

min(a$HEIGHT,na.rm=T)                         ## Min

## [1] 51.3

변수 만들기

1

HEIGHT, WEIGHT 변수를 이용해서 BMI 변수를 만들어보자.

a$BMI <- a$WEIGHT/((a$HEIGHT/100)^2)
a$BMI

##  [1]  29.18599  23.97353        NA  19.08418  22.53995  27.05627  22.75956
##  [8]  22.72571  23.40672  22.32128        NA  25.27620  22.33710  19.48452
## [15]  23.63281  45.48324  34.88653  27.16967  18.36706  16.09042  22.35622
## [22]  20.16475  22.31107  21.05703 189.23201  21.54913  22.65067  16.23245
## [29]  22.75374  35.71304        NA  12.69589  24.13891  19.75689  20.04338
## [36]  19.21222  20.67644  29.01123  23.17902        NA  21.73339  16.15029
## [43]  21.01780  23.50206        NA  23.80408  20.75674  24.75186  28.06813
## [50]  17.61528

2

25를 기준으로 BMI 분류 변수를 만들어보자

a$BMI_cat=0                           ## 0으로 된 변수만들고
a$BMI_cat[a$BMI>=25]=1                ## 25이상인 것만 1로 바꾼다.
table(a$BMI_cat)

## 
##  0  1 
## 40 10

a$BMI_cat=a$BMI>=25                  ## 방법 2: TRUE of FALSE 
table(a$BMI_cat)

## 
## FALSE  TRUE 
##    35    10

a$BMI_cat=ifelse(a$BMI>=25,"1","0")  ## 방법 3: ifelse문 이용
table(a$BMI_cat)

## 
##  0  1 
## 35 10

3

ifelse의 응용

a$BMI_cat=ifelse(a$BMI>=25,"overweight",ifelse(a$BMI<18.5,"underweight","normal"))
table(a$BMI_cat)

## 
##      normal  overweight underweight 
##          29          10           6

a

연속형, 범주형 변수 지정하기

1

a 자료의 대략적인 요약을 살펴보자.

summary(a)

##      CASEID           AGE            SEX                HEIGHT     
##  Min.   : 1.00   Min.   :15.00   Length:50          Min.   : 51.3  
##  1st Qu.:13.25   1st Qu.:27.50   Class :character   1st Qu.:158.4  
##  Median :25.50   Median :46.50   Mode  :character   Median :163.9  
##  Mean   :25.50   Mean   :44.32                      Mean   :162.6  
##  3rd Qu.:37.75   3rd Qu.:59.50                      3rd Qu.:172.5  
##  Max.   :50.00   Max.   :75.00                      Max.   :183.0  
##                                                     NA's   :3      
##      WEIGHT             DM            HTN             PH          TREATMENT  
##  Min.   : 40.00   Min.   :0.00   Min.   :0.00   Min.   :6.850   Min.   :0.0  
##  1st Qu.: 54.75   1st Qu.:0.00   1st Qu.:0.00   1st Qu.:7.122   1st Qu.:0.0  
##  Median : 60.00   Median :1.00   Median :0.00   Median :7.180   Median :0.5  
##  Mean   : 62.78   Mean   :0.66   Mean   :0.46   Mean   :7.176   Mean   :0.5  
##  3rd Qu.: 67.58   3rd Qu.:1.00   3rd Qu.:1.00   3rd Qu.:7.240   3rd Qu.:1.0  
##  Max.   :120.40   Max.   :1.00   Max.   :1.00   Max.   :7.370   Max.   :1.0  
##  NA's   :4                                                                   
##       BMI           BMI_cat         
##  Min.   : 12.70   Length:50         
##  1st Qu.: 20.16   Class :character  
##  Median : 22.54   Mode  :character  
##  Mean   : 26.80                     
##  3rd Qu.: 24.14                     
##  Max.   :189.23                     
##  NA's   :5

이상한 점은 변수의 class가 올바르게 설정되지 않았기 때문이다.

sapply(a,class)

##      CASEID         AGE         SEX      HEIGHT      WEIGHT          DM 
##   "integer"   "integer" "character"   "numeric"   "numeric"   "integer" 
##         HTN          PH   TREATMENT         BMI     BMI_cat 
##   "integer"   "numeric"   "integer"   "numeric" "character"

2

CASEID는 character로, SEX, DM, HTN, TREATMENT, BMI_cat 변수는 범주형 변수로 바꾸고자 한다.

a$CASEID <- as.character(a$CASEID)

cat_vars <- c("SEX","DM","HTN","TREATMENT","BMI_cat")

for(vn in cat_vars){
  a[,vn] <- as.factor(a[,vn])
}

sapply(a,class)

##      CASEID         AGE         SEX      HEIGHT      WEIGHT          DM 
## "character"   "integer"    "factor"   "numeric"   "numeric"    "factor" 
##         HTN          PH   TREATMENT         BMI     BMI_cat 
##    "factor"   "numeric"    "factor"   "numeric"    "factor"

summary(a)

##     CASEID               AGE          SEX         HEIGHT          WEIGHT      
##  Length:50          Min.   :15.00   F   :19   Min.   : 51.3   Min.   : 40.00  
##  Class :character   1st Qu.:27.50   M   :29   1st Qu.:158.4   1st Qu.: 54.75  
##  Mode  :character   Median :46.50   NA's: 2   Median :163.9   Median : 60.00  
##                     Mean   :44.32             Mean   :162.6   Mean   : 62.78  
##                     3rd Qu.:59.50             3rd Qu.:172.5   3rd Qu.: 67.58  
##                     Max.   :75.00             Max.   :183.0   Max.   :120.40  
##                                               NA's   :3       NA's   :4       
##  DM     HTN          PH        TREATMENT      BMI                BMI_cat  
##  0:17   0:27   Min.   :6.850   0:25      Min.   : 12.70   normal     :29  
##  1:33   1:23   1st Qu.:7.122   1:25      1st Qu.: 20.16   overweight :10  
##                Median :7.180             Median : 22.54   underweight: 6  
##                Mean   :7.176             Mean   : 26.80   NA's       : 5  
##                3rd Qu.:7.240             3rd Qu.: 24.14                   
##                Max.   :7.370             Max.   :189.23                   
##                                          NA's   :5

Subset 만들기

특정 조건을 만족하는 부분 데이터셋을 만드는 것은 지금까지 배웠던 것을 응용할 수도 있고 subset이라는 함수를 쓸 수도 있다.

a[a$AGE>=60 & a$SEX=='F',]                                    ## 행에 조건 써준다. 

subset(a,a$AGE>=60 & a$SEX=='F')                              ## 문법이 간단하여 이것을 주로 쓴다.

subset(a,a$AGE>=60 & a$SEX=='F', select=c("AGE","SEX","DM"))  ## select로 원하는 변수만 추려낸다.

Table과 Plot 예시

Table 1

대부분의 의학 논문에는 데이터의 특성을 나타내는 Table1이 포함되어야 한다.
Table1을 만드는 다양한 함수 중 tableone 패키지의 CreateTableOne 함수를 사용해보겠다.

# install.packages("tableone")
library(tableone)

tb1=CreateTableOne(vars = names(a[,2:9]),       ## 테이블에 들어갈 변수들
                   strata = "BMI_cat",          ## 그룹 변수(여러 개 가능)
                   factorVars = cat_vars,       ## 범주형 변수들
                   data = a,
                   includeNA = T                ## 범주형 변수에서 NA를 하나의 범주로 포함할 것인가?
                  )
tb1

##                     Stratified by BMI_cat
##                      normal         overweight     underweight    p      test
##   n                      29             10              6                    
##   AGE (mean (SD))     46.17 (15.78)  43.10 (22.98)  42.00 (21.12)  0.823     
##   SEX (%)                                                          0.805     
##      F                    9 (31.0)       5 (50.0)       2 (33.3)             
##      M                   19 (65.5)       5 (50.0)       4 (66.7)             
##      NA                   1 ( 3.4)       0 ( 0.0)       0 ( 0.0)             
##   HEIGHT (mean (SD)) 165.07 (8.25)  148.58 (34.76) 175.68 (2.45)   0.010     
##   WEIGHT (mean (SD))  59.98 (6.88)   76.81 (20.08)  49.96 (6.00)  <0.001     
##   DM = 1 (%)             20 (69.0)       5 (50.0)       4 (66.7)   0.554     
##   HTN = 1 (%)            17 (58.6)       3 (30.0)       2 (33.3)   0.211     
##   PH (mean (SD))       7.17 (0.10)    7.21 (0.07)    7.17 (0.18)   0.651     
##   TREATMENT = 1 (%)      15 (51.7)       4 (40.0)       2 (33.3)   0.636

Plots

R의 큰 장점 중 하나는 수려한 그래픽이다. 간단히 예제 데이터로 두 가지를 그려보자.
오늘은 데이터의 형태에 따라 상상하는 거의 모든 그래프를 그릴 수 있다는 것 정도만 기억하자.

1

library(ggplot2)

a<-na.omit(a)                                    ## 결측치가 없는 행만 추린다

ggplot(a, aes(x=HEIGHT, y=BMI, color=SEX)) +
  geom_point() +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE)

## `geom_smooth()` using formula 'y ~ x'

2

a$BMI_cat <- factor(a$BMI_cat, levels = c("underweight", "normal", "overweight"))

ggplot(a, aes(x=BMI_cat, y=PH, color=BMI_cat)) +
  geom_boxplot() +
  geom_jitter(color="black", size=0.4, alpha=0.9)

More examples

다양한 패키지, 함수, 옵션들을 활용하여 이런 표와 그림들을 그릴 수 있다.
그림을 잘 그리기 위해서는 데이터의 적절한 가공, 적절한 함수 선택, 가시성을 높여주는 적절한 옵션의 추가가 중요하다.

1

datatable 함수에 formatStyle로 디자인했다. Figure 3. Differences in dose distribution values
The values represent the differences between the values obtained by the techniques listed in the columns and those obtained by the techniques listed in the rows. Red indicates statistically significant-positive differences and blue indicates statistically significant-negative differences

2

ggplot, geom_line 함수를 사용했다. Figure 2. Dose-volume histograms (DVH)
(A) PTV, (B) IMN, (C) heart, (D) ipsilateral lung, (E) LAD, and (F) ipsilateral breast skin. Pale lines represent the DVH per patient per technique; thick lines represent the mean DVH curve of each technique
IMN, internal mammary lymph nodes; LAD, left anterior descending artery; PTV, planning target volume

3

ggboxplot 함수를 사용했다. Figure 4. Ipsilateral lung volumes receiving over 5, 20, and 30 Gy per patient per technique

4

ggplot, geom_bar, geom_errorbar 함수를 사용했다. Figure 5. Doses to organs-at-risk per technique, aimed at decreasing the mean heart dose

수고하셨습니다!
http://147.47.68.165:1111/jihyunbaek/lithium/
http://147.47.68.165:1111/doctorssi/ranitidine/